/ * METHOD FOR THE HIGHER-ORDER BLIND IDENTIFICATION OF 

MIXTURES OF SOURCES 

BACKGROUND OF THE INVENTION 
1. Field of the Invention 

The invention relates especially to a method of fourth-order and 
5 higher order self-learned (or blind) separation of sources from reception 
with N receivers or sensors (N>2), this method exploiting no a priori 
information on sources or wavefronts, and the sources being P 
cyclostationary (deterministic or stochastic, analog or digital sources with 
linear or non-linear modulation) and statistically independent sources, 
10 It can be applied for example in the field of radio 

communications, space telecommunications or passive listening to these 
links in frequencies ranging for example from VLF to EHF. 

It can also be applied in fields such as astronomy, biomedicine, 
radar, speech processing, etc. 
15 2. Description of the Prior Art 

The blind separation of sources and, more particularly, 
independent component analysis (ICA) is currently arousing much interest. 
Indeed, it can be used in many applications such as telecommunications, 
speech processing or again biomedicine. 
20 For example, in antenna processing, if signals sent from a certain 

number of sources are received at an array of receivers and if, for each 
source, the timing spread of the channels associated with a different 
receivers is negligible as compared with the timing symbol time, then an 
instantaneous mixture of the signals sent from the sources is observed on 
25 said receivers. 

The blind separation of sources is aimed especially at restoring 
the sources assumed to be statistically independent, and this is done solely 
on the basis of the observations received by the receivers. 

Depending on the application, it is possible to retrieve only the 
30 instantaneous mixture, namely the direction vectors of the sources. This is 
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1 the case for example with goniometry where said mixture carries all the 
information needed for the angular localization of the sources by itself: the 
term used then is "blind identification of mixtures". 

For other applications such as transmission, it is necessary to 
5 retrieve the signals sent from the sources: the expression used then is 
separation or again blind or self-learned extraction of sources. 

Certain prior art techniques seek to carry out a second-order 
decorrelation (as can be seen in factor analysis with principal component 
analysis (PCA). 

10 ICA, for its part, seeks to reduce the statistical dependence of the 

signals also at the higher orders. Consequently, ICA enables the blind 
identification of the instantaneous mixture and thereby the extraction of the 
signals sent from the sources, not more than one of which is assumed to be 
Gaussian. At present, this is possible only in complying with certain 

15 assumptions: the noisy mixture of the sources must be linear and 
furthermore over-determined (the number of sources P must be smaller 
than or equal to the number of receivers N). 

While Comon was the first to introduce the ICA concept and 
propose a solution, COM2 in the reference [1] (the different references are 

20 brought together at the end of the description) maximized a contrast based 
on fourth-order cumulants, Cardoso and Souloumiac [2], for their part 
developed a matrix approach, better known as JADE, and thus created the 
joint diagonalization algorithm. 

Some years later, Hyvarinen et al. presented the FastICA 

25 method, initially for real signals [3], and then in complex cases [4]. This 
method introduces a contrast-optimizing algorithm called the fixed-point 
algorithm. 

Comon has proposed a simple solution, COM1 [5], to contrast 
optimization presented in [6]. 
30 Although these methods perform very well under the 

assumptions stated here above, they may nevertheless be greatly disturbed 
by the presence of unknown noise, whether Gaussian or not, that is 
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1 spatially correlated and inherent in certain applications such as HF radio 
communications. 

Furthermore, as stated further above, the above methods are 
designed only to process over-determined mixtures of sources. Now in 
5 practice, for example in radio communications, it is not rare to have 
reception from more sources than receivers, especially if the reception 
bandwidth is great. We then have what are called under-determined 
mixtures (P> JV). 

Several algorithms have been developed already in order to 
10 process mixtures of this type. Some of them tackle the difficult problem of 
the extraction of sources [7-8] when the mixture is no longer linearly 
inverted, while others deal with the indication of the mixture matrix [7] [9- 
12]. 

The methods proposed in [9-11] exploit only fourth-order 

15 statistics while the method presented in [12] relies on the use of the 
characteristic second function of the observations, in other words on the 
use of non-zero cumulants of any order. As for the method used in [7], it 
relies on the conditional maximization of probability, in this case that of 
data conditional on the mixture matrix. 

20 While these methods perform well, they have drawbacks in the 

operational context. 

Thus, the method [9] is difficult to implement and does not 
ensure the identification of the direction vectors of sources of the same 
kurtosis. The methods [10] and [11] for their part cannot be used to 

25 identify the direction vectors of circular sources. The method [10], called 
S3C2, for its part confines the user in a configuration of three sources and 
two receivers, ruling out any other scenario. The method [7] authorizes the 
identification of four speech signals with only two receivers. However the 
samples observed must be temporally independent and each source must 

30 have a sparse density of probability. Finally, the method [12] is applicable 
only in the case of real sources, which is highly restrictive especially in 
digital communications. Furthermore, the algorithm depends greatly on the 
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•number of sources, and there is nothing today to prove that a poor 
estimation of this parameter will not impair the performance of the method. 
SUMMARY OF THE INVENTION 

The present invention offers a novel approach relying especially 
5 on the exploitation of the totality or practically the totality of the 
information proper to the direction vectors a p of the sources, contained 
redundantly in the matrix representing m = 2q order circular statistics of 
the vector of the complex envelopes of the signals at output of the receivers. 

The invention relates to a method for the blind identification of 

10 sources within a system comprising P sources and N receivers. It comprises 
at least one step for the identification of the matrix of the direction vectors 
of the sources from the information proper to the direction vectors a p of the 
sources contained redundantly in the m=2q order circular statistics of the 
vector of the observations received by the N receivers. 

15 The m = 2q order circular statistics are expressed for example 

according to a full-rank diagonal matrix of the autocumulants of the 
sources and a matrix representing the juxtaposition of the direction vectors 
of the sources as follows: 

0,* = Aq £m,s AqH (11) 

20 where £m,s = diag([C,* , y^ ,.m,C^J) is the full-rank diagonal matrix of the 
m-2q order autocumulants C p p p p ^'f p s des sources, sized (Px Pj y and where 

Aq = [a^ 9 ° (8) a* a p 0{q ° ® a p ] , sized (N<* x Pj and assumed to be of 

full rank, represents the juxtaposition of the P column vectors 

[a™® a']. 

L p p J 

25 The method of the invention is used in a communications network and/or 
for goniometry using identified direction vectors. 

The invention has especially the following advantages: 

• It enables the blind identification of instantaneous mixtures, both 
over- determined (where the number of sources is smaller than or 
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equal to the number of receivers) and under-determined (where the 
number of sources is greater the number of receivers) as well as the 
blind extraction of the sources in the over-determined case; 

• At the m=2q order, which is an even-parity value where q>2, the 
5 procedure called BIOME (Blind Identification of Over and 

underdetermined Mixtures of sources) can process up to P = Nta- 1 ) 
sources using the array with N different receivers, once the m order 
autocumulants of the sources have the same sign; 

• An application of the fourth order method known as ICAR 
10 (Independent Component Analysis using Redundancies in the 

quadricovariance), enables the blind identification of over-determined 
(P<N) instantaneous mixtures of sources and their blind extraction, 
in a manner that proves to be robust in the presence of a spatially 
correlated unknown Gaussian noise, once the sources have same- 
15 sign kurtosis (fourth-order standardized autocumulants ); 

• an application of the BIOME method at the sixth-order level, called 
BIRTH (Blind Identification of mixtures of sources using 
Redundancies in the daTa Hexacovariance matrix), enables the blind 
identification of instantaneous mixtures, both over-determined {P<N) 

20 and under- determined (P > N), of sources, as well as the blind 

extraction of the sources in the over-detern^ned case. The BIRTH 
method has the capacity to process up to N 2 sources from an array 
with N different receivers, once the sixth-order autocumulants of the 
sources have a same sign. 

25 BRIEF DESCRIPTION OF THE DRAWINGS 

Other characteristics of advantages of the method according to the 
invention shall appear more clearly from the following description from a 
non-restrictive example of an embodiment and the appended figures, of 
which: 
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• Figure 1 is a drawing exemplifying an m-order implementation of the 
method, 

• Figures 2, 3 and 4 show results of simulation of a fourth-order 
implementation of the method, 

5 • Figures 5, 6 and 8 show results of simulation of a sixth-order 
implementation of the method. 

The following examples are given for the identification and /or 
extraction of sources in an array comprising an array antenna comprising N 
receivers. It is assumed that this antenna receives a noisy mixture of 
10 signals from P statistically independent sources for example in narrow band 
(NB). 

On the basis of these assumptions, the vector x(t\ of the complex 
envelopes of the signals at output of the receivers is written, at the instant t 

p 

*M-£ s P (tia P A 44 + i<<) (1) 

P =\ 

15 • where v(t) is the noise vector, assumed to be centered, Gaussian, 

spatially white and unknown, 
• s p (t) and a p correspond respectively to the complex, narrow-band, 

cyclostationary and cycloergodic envelope BE (with a possible residue of 

a carrier as the case may be) and to the direction vector of the source p, 
20 • s(t) is the vector whose components are the signals s p (t) and A is the 

matrix (N x P) whose columns are the vectors a p . 

Furthermore, the method generally described here below for the 
m = 2q (q > 2) order uses the following assumptions, numbered Hi_4: 

HI : At any instant t, the sources with complex values s p (t) are 
25 cyclostationary, cycloergodic and mutually decorrelated at the m order; 

H2 : At any instant t, the components v n (t) of the noise are stationary, 
ergodic, Gaussian and circular; 
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* H3 : A any instant t, s(t) and v(tj are statistically independent; 

H4 : the m order autocumulants of the sources are not zero and have the 
same sign. 

With the above assumptions, for a given even-parity order m = 
5 2q, the problem of the blind identification of instantaneous mixtures of 
sources consists in finding the matrix A through the exploitation of certain 
m order statistics of the observations. This matrix A is found to the nearest 
trivial matrix (a trivial matrix has the form A IT where A is an invertible 
matrix and II is a permutation matrix). 

10 The blind separation (extraction) of sources consists especially in 

determining the separator that is linear and invariant in time (LIT), W, with 
a dimension (Nx F), the output vector of which has a dimension (Px 1), 

iM- WHx(t) (2) 

corresponds, plus or minus a trivial matrix, to the best estimate, S(/) , of 
15 the vector s(t), where the symbol H signifies a conjugate transpose. 

The separator W is defined to the nearest trivial matrix inasmuch 
as neither the value of the output power values nor the order in which they 
are stored changes the quality of restitution of the sources. 

Before any explanation of the steps of the method according to 
20 the invention, a few reminders are given on the statistics of the 
observations. 

Statistics of the observations 

The method according to the invention uses especially 2, 4, 6, 
m even-parity circular statistics of the observations. 
25 According to the prior art described in the reference [15], the 

expression of the m order cumulants as a function of the moments of an 
order lower than m can be simplified. 

Let Gfc"£ x be a scalar quantity with a complex value depending in the q 
lower indices d, e and the q higher indices /, g having values in {1, 2, 
30 N\. The quantity G^f x then verifies the following three symmetries: 
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• any permutation between the lower indices of G^^f does not modify 
the value of G£"'J X : for example, for q = 2, G{ d 8 x = Gf* , 

• any permutation between the higher indices of G^ 'f x does not 
modify the value of Gj^f s : for example, for q = 2, G^{ x = Gjf* , 

• permutating all the higher indices with all the lower indices of G^;f x 

has the effect of conjugating the value of G f £ * x : Gf '^ g x = (g^^)* . 
Furthermore, the following notation is adopted: the quantity 
i r ]Gj^f x Gfc' f x ... G"'[^ x designates the linear combination of the r possible 
and distinct products (modulo the three symmetries described here above) 
of the type G^J X G J h ^ x ...G";^ xi weighted by the value 1. Each of the r 

products is built from the product G^'"f x G£"J* X •••G / ";;;j£ in using and 
combining the following two rules of permutation: 

• a lower index of one of the terms of the product G^'f x Gfc'f x ... G? 7 ""° x 
permutates with a lower index of another term of the same product 
GjZll x ^^x---^t'm,x to S ive another (distinct) product: for example, 
for (gi, q 2 ) = (2, 2), G^f x G J h j x gives as other distinct products (modulo 
the three symmetries described here above) G{f x G J d j x , G$f x Gj; k ex 

and G f * 8 G J ' k 
ana u h j x v d e x , 

• a higher index of one of the terms of the product 
Gj^£ x Gfc''j x ... G"; '* x permutates with a higher index of another term 

of the same product G£;;£ x Gfcj x ... G"'^ x to give another (distinct) 
product: for example, for (qi, 92) = (2, 2), G^f x G J h f x gives the following 
as other distinct products (modulo the three symmetries described 
here above) Gj;f >JC G{f x , G% x G& and G>*„ G{« x , 

in order to obtain the totality of the r possible and distinct products (modulo 
the three symmetries described here above). 



9 

' The following example, where (q h q 2 ) = (2, 1) et r = 9, illustrates this 
notation: 

[9] Gf h eiX G' Lx = G&G^+G&Gi, + Gf* x G' dx + G£, G*. 
+ % + G* x + Gj'j x Gl + G*i x Gl + G^ x Gl x (3) 

5 m order statistics 

In the case of potentially non-centered, stationary or 
cyclostationary sources, the m order circular statistics of the vector x(tj, 
given by (1), can be written: 

0 

where the symbol < /(/) > f = \\m{\IT)\ T J r ) 2 f(t)dt corresponds to the 

7"->oo 

operation of time-domain averaging, in continuous time, of f(t) on an infinite 
horizon of observation. 



15 C^t 1 '^^ (0,...,*, (/),*, (0*,*, (0*.-,*,- (0*) (5) 

l Y l V~ J q* \ l 2 1 <1 l q+\ l q+2 l ™ 

where q terms are conjugate and q terms are non-conjugate. 

The m order statistics described by the expression (5) are said to be 
circular because the m order cumulant Cum{xd(f), xw(*),..., Xg(t)* y 
Xh( *)*,..., Xi(t)*} is computed from as many conjugate terms (Xg(t)*, Xh(t)*, 
20 x,(f)*) as non-conjugate terms (xd(*), Xe(tj, xy(t)). 

The m order circular cumulants may be expressed as a function 
of the lower- than- m order moments as follows. 



25 



10 

Let M^ Vtq f r " Jm (t) and C^f'^f 1 '"^ \t) be the m order 

moments and circular cumulants associated with the observation vector 
x(t), defined by : 

M >+lV2^% ) = E[ x (t) je (ft x (t)>JC (t) * )JC (t) * 5 .. (t) *j 

5 

G 

/Vir/* s ={s* s (l)U5* I (l)U.-.USi tX (*)} designates the u g v th partition, among 
a Gk possible partitions, of « k » subsets Sp x Juty 9 "'' 1 ™ w here 1 < r * s * ... 

* t < q et q+1 < u * v * ... * w < m, such that Part% x ■ v . As for the 
1 0 union operator, written as U , it verifies 

In practical terms, two great classes of estimators may be used to 
estimate the above statistics: in the case of stationary sources, it is possible 
to use a non-skewed and consistent empirical estimator for potentially 

15 centered, ergodic, stationary sources whereas, in the case of potentially 
non-centered cyclostationary sources, it is necessary to use what is called 
an exhaustive, non-skewed and consistent estimator for potentially non- 
centered, cycloergodic and cyclostationary sources. This exhaustive 
estimator is, for example, determined according to an approach described in 

20 the references [13-14]. 

Arrangement and storage of m order statistics 

As shown here above, the statistics , for 1 < d, e f, 

tf,e,...,/,x' ' Jt 

g f h, k < N are functions with m = 2q inputs (m being an even-parity 
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. value), it is then possible to arrange them in a matrix (N* x Nv) that will be 
named 

Explicitly, the quantity C^'j x is located at the ith row and at the jth 

column of the matrix Cm,* in writing i = N[...N\N(d-l)+e-l]+ ...]+/c t j = N[ 
5 ...N[N[g-l) + ...]+/ 

Arrangement and fourth-order statistics 

In the case of centered stationary sources, the fourth-order 
circular statistics of the vector x(*), given by (1), are written as follows: 

Cj£ = Cum{^(i), xJM, xAV> - " M dM M[* - [2]M^ M* x (4)' 

10 In practical terms, these statistics may be estimated by using a 

non-skewed and consistent empirical estimator for centered, ergodic, 
stationary sources. 

In the case of potentially non-centered cyclostationary sources, 
the fourth-order circular statistics of the vector to be taken into account 
15 are written as follows: 

CJ£ = < Cumte(t), xe(t), xAt)\ xJW\ > c = < M^ x (t) > c - < [*]M dx (t)M f J(t) 

>c~< M dfC Jt)M^(t) > c -< [2)Mi % (t)Ml(t) > c + 2 < M„, x (f)M c%x (t)M^(t) 

> c + 2 < A/^WAffWA/^tl) > c + 2 < [4]M^(t)M/ (t)M^(t) > c - 6 < 

M J Jt)M e JQM x f (t)M*(1)> c (5)' 

20 In practical terms, these fourth-order statistics may be estimated 

by using the estimator known as the exhaustive, non-skewed and 
consistent estimator for cyclostationary, cycloergodic and potentially non- 
centered sources. This exhaustive estimator is described in 13-14]. 

As shown here above, the statistics C d f x , for 1 < d, e, f } g < N are 

25 four-input functions. It is possible then to arrange them in a (N 2 x N 2 ) 
matrix that will be called a quadricovariance matrix Q x given for example in 
the reference [13-14], 
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' Arrangement and sixth-order statistics 

In the case of centered stationary sources, the sixth-order 
circular statistics of the vector x(£), given by (1), are written as follows: 

CJ# iB - CumWt), xJto, xAt), xg(1)* t XhW, *{<)*} = M^ fx - [3]Ml f x M h / - 
5 [9)M% X M' Lx - [3)M dex Mff + 2[9]M^ X A/^ + 2[6]A/* X A/* s A^ (6) 

In practical terms, these statistics may be estimated by using a 
non-skewed and consistent empirical estimator for centered, ergodic, 
stationary sources. 

In the case of centered cyclostationary sources, the sixth-order 
10 circular statistics of the vector x(tj given by (1), are written as follows: 

C#k = < CumWfl, xjfa xM* 9 xhW, x(fl*} > c = < M^Jtf > c - < 
[3] (1) Af ^ (t) >c - < [9] M^ x (t) M* f x (f) > c - < [3] (t) (£) > c + 2 < 

[9] Af ^ (<) M* , (£) A/," (t) > c + 2 < [6] A/* M (I) M h c x (t) M' f x (t) > c (7) 

In practical terms, these sixth-order statistics may be 
15 estimated by using an estimator called an exhaustive, non-skewed and 
consistent estimator for cyclostationary, cycloergodic and centered sources.. 

As shown here above, the statistics C£*'* x , for 1 < d, e, f, g, h, k < 

N are six-input functions. It is possible then to arrange them in a (iV 3 x N 3 ) 
matrix that will be called a hexacovariance matrix H x . 

20 

Principle implemented in the method according to the invention 

The method according to the invention uses especially a property 
of multilinearity of the cumulants and the Gaussian nature of noise which 
take the form of the following matrix expression 

25 Cm,* = i) ® A*] Cm.s [A%- D ® A*}" (10) 

where Cm,* and Cm,s are the matrices of the m order statistics defined earlier, 
having respective sizes x and (Pu P<?), and being associated with 
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' the vectors x(t) and s(t) where <A®(<? - *> corresponds to an adopted notation 
defined thus: the matrix designates the matrix B raised to the power (in 
the sense of the Kronecker product) fc, i.e. in taking the Kronecker product 
Ef*k = B®Bg>. ..®B , in writing = 1. 

k times 

5 The Kronecker product may be recalled here: let A and B be two matrices 
respectively sized (L A x C A ) and (L B x Cb). The Kronecker product D= A® B 
is a matrix sized (L A L B x Ca C b ) defined by D = (Aij B) l<i< LA, l<j< CA. 

Without departing from the framework of the invention, other 
modes of expression associated with other modes of arrangement of the 
10 cumulants may be used: 

Cm,* = [A®*] Cm,s,l[A®<l} » (10a) 

where 1 is chosen arbitrarily between 1 and q and where Cm,s, r. is the matrix 
of the m = 2q order statistics of s(t) associated with the index 1 chosen. 
Each expression conditions the number of sources potentially identifiable 
15 from a given array. 

Here below in the description, the analysis is given in using the expression 
of the relationship (10). 

Inasmuch as the sources are independent, the matrix of the m 
order statistics associated with the sources, C m $i is a diagonal matrix. 
20 However, it turns out to be not a full-rank matrix. The method according to 
the invention considers a matrix determined from a full-rank diagonal 
matrix of the autocumulants and from the matrix representing the 
juxtaposition of the P column vectors relative to the direction vectors of the 
sources: 

25 Cm,x = A q fasA c p (11) 

where £m, s " diag{[C^^ ,...,C^J) is the full-rank diagonal matrix of the 
m = 2q order autocumulants C^ ''^ s from the P sources, sized (Px F), and 
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fl, 09 a, a p 09 J sized (N<* x F) and assumed to 

be of full rank, represents the juxtaposition of the P column vectors 

[CL P ® fa 1 ) ® Op* ]. Furthermore, we assume that the matrix 

\n^ q ~ 2) ten* n ^ q ~ 2) (tin *1 
Ah = l a \ —- a p ^ a p J, sized (NW x P\, is also a full- 

5 rank matrix. 

The method according to the invention enables the 
advantageous exploitation and extraction, for example, of the totality of the 
information proper to the direction vectors a p of the sources, redundantly 
contained in the matrix of the m = 2q order circular statistics of the 
10 observations vector x(t), Cm,* and more particularly in the matrix Aq. 

The method comprises, for example, the steps described here 
below. The samples of the vector x(f) are assumed to have been observed 
and the matrix Cm,x is assumed to have been estimated from these samples. 

Step 1: Singular value decomposition of the matrix C m jr 

15 This step computes the square root C mx f/2 of the full-rank matrix 

Cm,x, for example through the eigenvalue decomposition of the Hermitian 
matrix C m>x - E s L s E S H where L s and E s are respectively the diagonal matrix 
of the P greatest (in terms of absolute value) real eigenvalues of C m , x and the 
matrix of the associated orthonormal eigenvectors. This step shows the 

20 relationship existing between O m J /2 and Aq : 

C2 = E\L\ n = A C U2 V H = [a^- x) ®a\ ® a ^ VI V H (12) 

m,x s | s | q m,s LI I p p J? m,\ » ' 

where V is a unit matrix, sized (P x P\> unique for L s and E s as given 
matrices, and where J / 2 , £*ii 2 are square roots respectively of \L S \ and 
£n,s ( | . | designates the absolute value operator). 

25 L s and E s are, for example, respectively the diagonal matrix of the P greatest 
(in terms of absolute value) real eigenvalues of Cm, x and the matrix of the 
associated orthonormal eigenvectors. 
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For a full-rank matrix Aq, it is possible to ascertain that the 
hypothesis (H4) is equivalent to assuming that the diagonal elements of L s 
are non-zero and have the same sign. 

Step 2 

5 This step consists of the extraction, from the matrix C l m '] = [ r i T , 

IV] T , of the N matrix blocks T n : each block T n sized (Nta- l ) x F) is 
constituted by the JVta- i) successive rows of C l m n x starting from the a N (<?- *) 

(n-l)+l"th row. 
Step 3 

10 This step entails the building of the N(N-l) matrices 0 n i,n2 

defined, for all 1 < rt\ * n 2 < N, by © n i,n2 = T * } T n2 where # designates the 
pseudo-inversion operator. 

In noting for all 1 < n < N, O n = diag([A n i, Ara, .... Anp]) where Aij is the 
component of A located on the ith row and jth column, there is equality T n = 

15 A iq _ }) ® n £ l J*V H for all 1 < n < N, and the fact that the matrix V jointly 

diagonalizes the N(N-1) matrices ® n x t n2 - T *T = VC' U2 <&~ 1 ® r U2 V H , 
which, it may be recalled, is sized (Px Pj. 
Step 4 

This step consists in determining the matrix V so i, resolving the 
20 problem of the joint diagonalization of the N(N-l) matrices 0ni,n2 for example 
in using a method of diagonalization described in the reference [2]. The 

matrix C™\ Vsoi, where V so i = V Tis a unitary matrix jointly diagonalizing the 

matrices 0ni,n2 to the nearest unitary trivial matrix T , is an estimate of the 
matrix A? to the nearest trivial matrix. 

25 Different methods, known to those skilled in the art, enable the 

extraction from Aq of an estimate A; A of the mixture matrix A, 

Step 5 
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Step 5A 

One procedure consists, for example, in taking the average of the 
K = ffiq- i) blocks (Z fc )* sized (N x P) (for all 1 < k < Nto- D (the block I* is 
constituted by N successive rows of Aq = [Ei T , I,K T ] r starting from the U N 
5 (fc-lJ+Fth row), or else in starting only one, for example (Zi)*. This approach 
enables the estimation, in any order and excepting an amplitude, of the P 
direction vectors a p and therefore the mixture A matrix to the nearest trivial 
matrix. 

Step 5B 

10 Another step consists, for example, for each of the P columns b p 

oi^ia^^a; a/^fca/Jin, 

• extracting the K = M<t~ 2 ) vectors b P (k) stacked one beneath the other such 
that: 

bp = [a* (q ~ X) (8) fl/]= [b P (l)T b p (2)T b p (K)T]T (i 4) then 

15 • converting said column vectors b p (k) = (Aip... Aj p ) \ci p (8> a ] sized (N 2 

x 1) into a matrix B p (k) = (Ap... Aj p ) [a p afi] (where 1 < i, j < N) sized (Nx N) 
and 

• jointly decomposing these K = M^ 2 ) matrices into singular values 
(singular value decomposition or SVD) : the eigenvector common to the K 

20 matrices B p (k) and associated with the eigenvalue that is the greatest (in 
terms of modulus) is therefore a column vector of the matrix A. It must 
be noted that the quantity (Aip... A jp ) is in the present case the product of 
(q - 2) components of A. 

This step of processing on the P columns b p of Aq, enables the 
25 estimation, in any order and plus or minus one phase, of the P direction 
vectors Op and therefore, to the nearest trivial matrix, the mixture matrix A. 

Step 6 

The mixture matrix A representing the direction vector of the 
sources contains, by itself, the information needed for the angular 
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localization of the sources. In this context, from the estimation of the 
different columns of A, it is possible to implement an arbitrary method of 
goniometry exploiting this information. Such a method is presented for 
example in the document [18]. 

5 Step 7 

To estimate the sources vector s(t) in an over- determined context 
(i.e. when P < N), the method applies an LIT type filter to the observations 
x(t) explicitly using the estimation of the mixture matrix A. It is possible, for 
example, to choose the FAS filter described in the reference [17], which is 
10 optimal in the presence of decorrelated sources. 

The method comprises, for example, a step 0 which consists of 

A 

the building, from the different observation vectors, x(tj, of an estimate C m 

of the matrix of statistics Cm,* of the observations, according to the method 
15 given earlier. In this case, the steps 1 to 6 of the method are implemented 

on the estimate C of the matrix. 

Criterion of performance 

According to one alternative embodiment, the method comprises 
20 a step using a normal criterion of performance for the evaluation of the 
blind identification of mixtures. This criterion is not global and enables the 
evaluation of the quality of identification of each direction vector estimated: 
it is then possible to compare two distinct methods of blind identification 
with respect to each direction vector, and hence to each source. This 
25 criterion is the extension, to the blind identification of mixtures, of the 
criterion based on SINR (signal-to-interference-plus-noise ratio) given in the 
reference [17] introduced for the blind extraction of sources. It is a P-uplet 
described by 

D(A, A; A ) = (<xi, a2, otp) (15) where 
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<*p = m«n [d(Op, d^)] (16) 

and where d{u,v) is the pseudo-distance between the vectors u and v, 
defined by 

d(«,H = l-|(«,v)| 2 |H|- 2 ||v||- 2 (17) 

5 It may be noted that < . , . > designates the scalar product defined for two 
vectors of a same dimension. 

The method described for the m = 2q order application can be 
applied especially to the fourth-order and sixth-order statistics, for example 
according to the examples given here below. 
10 Application of the method for the blind separation of fourth-order 
sources 

An alternative embodiment of the method known as ICAR (Independent 
Component Analysis using Redundancies in the quadricovariance) exploits 
the m = 4 (q = 2) order statistics, corresponding to the matrix of circular 
1 5 quadricovariance of the Cm,x, written as This method enables the blind 
identification of the instantaneous mixture A or the blind extraction of the 
sources s(tj when N > P, in other words only when the mixture is over- 
determined. 

The model (1) is assumed to be verified along with the fourth- 
20 order hypotheses Hm. 

The method called ICAR exploits especially the expression (11) which, when 
written for fourth-order statistics, is expressed as follows : 

Q X = A 2 £uA 2 h (18) 

where = diag([C 1 1 ; 1 1 t5 ,...,C££ S ]) is the full-rank matrix of the fourth-order 
25 autocumulants C£/? iS of the sources, sized (P x F) t and where Ai - 
[a ] ® a x a p <S>a p ] sized (N 2 x F) and assumed to be full-rank, 

represents the juxtaposition of the P column vectors [ Q ® a ] . 
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Furthermore, assuming that the mixture A matrix sized (N x F), is also a 
full-rank matrix. 

The method performs, for example, the steps 0 to 5 described in 
the case of the m=2q order application, in using the following parameters: 
5 Cm,x = Q x and =£m 

In this example of an implementation, the method may also 
include a step 0 which consists of: the building, from different observation 

vectors x(t), of an estimate Q of the matrix of quadricovariance Q x of the 
observations. The steps 1 to 6 are then carried out on this estimated value. 

10 Examples of results obtained by applying the method to fourth-order 
statistics 

Figures 2, 3 and 4 show a graph in which the x-axis corresponds 
to the number of samples L and the y-axis to the performance, the results 
of the simulation presenting the performance of the method according to 
1 5 the invention in a fourth-order application presented here above of ICAR (in 
implementing the step 5B), ICAR2 (in implementing the step 5A) and 
methods for the blind separation of sources (COM1, COM2, JADE, FastICA) 
known to those skilled in the art. The conditions of simulation are the 
following: 

20 • It is assumed that signals from P = 3 non-filtered sources, namely one 
BPSK source and two QPSK sources, are received on a circular array of 
N = 5 receivers such that R/X = 0.55 (with R and X the radius of the 
array and the wavelength) and such that the signal-to-noise ratio, SNR, 
is equal to 20 dB for each source. 

25 • The noise is Gaussian and spatially non-correlated. 

• The three sources are baseband sources and their time symbol is chosen 
to be equal to the sampling time. 

• The criterion used to obtain the best appreciation of the results of 
extraction from the source p for a given method is the maximum 

30 signal-to-interference-plus-noise ratio associated with the source p, 
better known as SINRM [17]. It may be compared with the optimum 
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SINRMp computed by using not the estimated mixture matrix but, on the 
contrary, the exact mixture matrix as well as the exact statistics of the 
observations. It is this comparison that is presented in figures 2-4. 

More particularly, figure 2 represents the SINRM of the source 1 
5 associated with ICAR, ICAR 2 along with the most efficient methods 
currently being used for the blind separation of sources such as JADE, 
COM1, COM2 and FastlCA. 

Figure 3 shows the SINRM of the source 2 for the same methods 
(ICAR, ICAR 2 , JADE, COM1, COM2, FastlCA and figure 4 shows those of the 
10 source 3. 

In each of the figures, it can be seen that the two methods ICAR 
and ICAR 2 perform very well and are slightly more efficient than the other 
methods JADE, COM1, COM2 and FastlCA. As for the FastlCA algorithm, 
its best performance relates to the source 2 for which it converges 

15 completely from 550 samples onwards. 

As for the difference in results between ICAR and ICAR 2 in this 
configuration it proves to be negligible as compared with the difference in 
performance between the ICAR methods and the JADE, COM1, COM2 and 
FastlCA methods which, however, are very good. 

20 Application of the method to sixth-order statistics 

According to another alternative embodiment, the method uses 
sixth-order statistics. This variant known as BIRTH (Blind Identification of 
mixtures of sources using Redundancies in the ddXa Hexacovariance matrix) 
enables the identification from an array of N receivers, of the direction 

25 vectors of at most P = N sources (in the case of an array with different 
receivers). It also enables the extraction of at most P = N sources for which 
the direction vector's are explicitly identified . 

The BIRTH method uses certain sixth-order statistics, stored in 
the hexacovariance matrix C6,x designated H x . Thus, this alternative 

30 implementation fully exploits the information proper to the instantaneous 
mixture A of the sources, contained in H x , especially through an artful 
writing of H x relative to the direction vectors of the sources, this being done 
by means of the property of multi-linearity of the cumulants: 
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Hx = Aa&sAa" (20) 

where &s = diag([Cl;l;l s ,...,Cp;p;p s ]) is the full-rank matrix of the sixth-order 
autocumulants C£#£$ of the sources, sized (Px F), and where 

A3 =[tf® 2 (8)^* a® 2 ® a*], sized (i\P x P) and assumed to be a full- 

5 rank matrix, represents the juxtaposition of the P column vectors 
r<3 02 ® a ] = la ® a ® a* • Furthermore, we assume that the matrix 

p P J L p p p J 

Ai =[a { ® a* a p ® a*] , sized (N 2 x F), is also a full-rank matrix. 

The sixth-order method comprises the steps 1 to 6 and the step 0 
described here above in using the following parameters: C ra ,x = H x et m=3. 

10 

Simulations 

Figures 5, 6, 7 and 8 show a graph in which the x-axis 
corresponds to the number of samples L and the y-axis corresponds to 
performance, namely the efficient operation of the method according to the 
15 invention in a sixth-order application. The following are the conditions for 
obtaining the curves: 

• Let it be assumed that P = 3 statistically independent sources, more 
particularly two QPSK sources and one BPSK source, all three being 
unfiltered, are received on a linear array of N = 2 receivers such that R / 

20 X = 0.55 (with R and X respectively being the radius of the array and the 
wavelength). 

• The three sources, assumed to be synchronized, have the same signal- 
to-noise ratio, written as SNR and being equal to 20 dB for each source 
with a symbol time that is four times the sampling time. 

25 • The BPSK source is chosen in baseband while the two QPSK sources 
have carriers respectively equal to half and one-third of the sampling 
frequency. 
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• The mixture matrix A is chosen so that the column vectors of the matrix 
A3 are linearly independent. The noise for its part is Gaussian and 
spatially non-correlated. 

• The instantaneous mixture of noisy sources is considered to be an over- 
5 determined mixture because the number of sources is greater than the 

number of receivers. The algorithms JADE, COM1, COM2, S3C2 known 
to those skilled in the art and the method of the invention in a sixth- 
order application are implemented for the blind identification of the over- 
determined mixture A. 
10 • The performance criterion a Pi defined by the equation (18), is computed 
on 200 operations, and this is done for each source p (1 < p < 3): it will 
thus enable the comparison of the five methods. 

According to the above assumptions, figure 5 shows the 
variations in the quantity ct3 resulting from the algorithms JADE, COM1, 
15 COM2, S3C2 and the method according to the invention, BIRTH, depending 
on the number of samples. The method according to the invention, unlike 
the prior art methods, makes it possible to identify the directional vector in 
question. 

Figure 6 gives a view, in the same context, of the variations of the 
20 triplet D(A, A\ A ) = (cti, cc2, (X3), associated with the method according to the 
invention in a sixth-order application as a function of the number of 
samples. The three coefficients a P rapidly decrease to zero as and when the 
number of samples increases. 

Figure 5 shows the variations in the quantity 013 resulting from 
25 the prior art methods JADE, COM1, COM2, S3C2 and the method 
according to the invention depending, this time, on the signal-to-noise ratio 
(SNR) of the source 3. The method BIRTH is fully successful in identifying 
the direction vector of the source 3 even for a low value of SNR. 



Finally, let it be assumed that the above P ~ 3 sources are received on a 
30 circular array of N = 3 receivers such that R/X = 0.55. 
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Figure 8 then shows the variations of the quantity 013 resulting 
from the algorithms JADE, COM1, COM2 and BIRTH has a function of the 
number of samples: the BIRTH algorithm works in an over-determined 
context, namely in the context where the number of sources is smaller than 
5 the number of receivers, and although sixth order cumulants of the 
observations must be estimated, the speed of convergence of BIRTH is the 
same magnitude is that of the methods referred to further above. 
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